'''
Created on Dec 2, 2010

@author: oabalbin
'''

def unified_genotyper_multiSample(ref_genome,list_bam_files, outfile, use_mem, num_cores, 
                           path_to_gatk, interva_list=None, call_parameters=None):
    '''
    This command will call variants using a list of bam files
    The bam files represent samples and are assumed to be properly
    processed
    java -jar -Xmx8g /exds/sw/bioinfo/gatk/GenomeAnalysisTK-1.0.4705/GenomeAnalysisTK.jar -T UnifiedGenotyper 
    -R /exds/projects/alignment_indexes/gatk/hg19/hg19.fa 
    -I /exds/users/oabalbin/projects/snps/exomes/aM18/test/s_3_12_sequence.hg19.aln2.rmdup.sorted.realigned.bam 
    -o /exds/users/oabalbin/projects/snps/exomes/aM18/test/s_3_12_sequence.snps.raw.vcf -stand_call_conf  [50.0] 
    -stand_emit_conf 10.0 -dcov [50] -L my.interval_list
    '''
    
    
    default_stand_call_conf=str(50.0)
    default_stand_emit_conf=str(10)
    default_dcov=str(50)
 
    # set up sample lists   
    input_files=[]
    for bamfile in list_bam_files:
        input_files += ['-I']+[bamfile]
    
    # call in specific sites
    if interva_list is not None:
        # call in specific sites
        call_on_this_regions=['-L',interva_list]
    else:
        #generate calls in all sites
        call_on_this_regions=['']

    # extra parameters
    if call_parameters is not None:
        stand_call_conf=str(call_parameters[0])
        stand_emit_conf=str(call_parameters[1])
        dcov=str(call_parameters[2])
    else:
        # Note 12-7-10. Check this parameters
        stand_call_conf=default_stand_call_conf
        stand_emit_conf=default_stand_emit_conf
        dcov=default_dcov

            
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'    
    args1 = ['java','-Xmx'+str(use_mem)+'m', '-jar',gatk_command,'-T','UnifiedGenotyper',
            '-l', 'INFO','-nt', str(num_cores),'-R', ref_genome, '-o', outfile]
    args2=  ['-stand_call_conf', stand_call_conf,'-stand_emit_conf',
             stand_emit_conf, '-dcov',dcov, '-A','AlleleBalance','-U', 'ALLOW_UNINDEXED_BAM']
    
    args = args1+input_files+args2+call_on_this_regions
    
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command


def unified_genotyper_pairedSample(ref_genome,list_bam_files_normal, list_bam_files_tumor,
                                   outfile, use_mem, num_cores, path_to_gatk, 
                                   interva_list=None, call_parameters=None):
    '''
    Important Note:
    In this mode all normal/tumor bam files are joined together
    and calls is made by comparing ALL reads for normal vs tumor
    Because of that, this mode should be used with for example
    a tissue cohort, e.g Prostate, but nor with mixed cohorts,
    e.g, breast plus Prostate cohorts data. 
    '''
    
    default_stand_call_conf=str(50.0)
    default_stand_emit_conf=str(10)
    default_dcov=str(50)
 
    # set up sample lists   
    normal_files=[]
    for bamfile in list_bam_files_normal:
        normal_files += ['-I:normal']+[bamfile]
    tumor_files=[]
    for bamfile in list_bam_files_tumor:
        tumor_files += ['-I:tumor']+[bamfile]
    
    # call in specific sites
    if interva_list is not None:
        # call in specific sites
        call_on_this_regions=['-L',interva_list]
    else:
        #generate calls in all sites
        call_on_this_regions=['']

    # extra parameters
    if call_parameters is not None:
        stand_call_conf=str(call_parameters[0])
        stand_emit_conf=str(call_parameters[1])
        dcov=str(call_parameters[2])
    else:
        # Note 12-7-10. Check this parameters
        stand_call_conf=default_stand_call_conf
        stand_emit_conf=default_stand_emit_conf
        dcov=default_dcov
    
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    args1 = ['java','-Xmx'+str(use_mem)+'m', '-jar',gatk_command,'-T','UnifiedGenotyper',
            '-l', 'INFO','-nt', str(num_cores),'-R', ref_genome, '-o', outfile]
    args2=  ['-stand_call_conf', stand_call_conf,'-stand_emit_conf',
             stand_emit_conf, '-dcov',dcov, '-A','AlleleBalance',
             '-U', 'ALLOW_UNINDEXED_BAM']
    
    args = args1 + normal_files + tumor_files + args2 + call_on_this_regions
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command




def indelGenotyper_singleSample(ref_genome,list_bam_files, outfile_vcf, outfile_bed, use_mem, num_cores, 
                                     path_to_gatk, interva_list=None, call_parameters=None):
    '''
    java -jar /path/to/GenomeAnalysisTK.jar \
     -T IndelGenotyperV2 \
     -l INFO \
     -R reference.fasta \
     -I sequencing.data.bam \
     -bed my.brief.output.bed        \
     -verbose my.detailed.output.txt \
     -o my.output.vcf \
     --refseq /path/to/refseq.rod \
     -L chr1
    '''
    if interva_list is not None:
        # call in specific sites
        call_on_this_regions=['-L',interva_list]
    else:
        #generate calls in all sites
        call_on_this_regions=['']
    
    if call_parameters is not None:
        stand_call_conf=call_parameters[0]
        stand_emit_conf=call_parameters[1]
        dcov=call_parameters[2]
    else:
        # Note 12-7-10. Check this parameters
        stand_call_conf=50.0
        stand_emit_conf=10
        dcov=50


    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    args1 = ['java','-Xmx'+str(use_mem)+'m', '-jar',gatk_command,'-T','IndelGenotyperV2',
            '-l', 'INFO','-R', ref_genome, '-o', outfile_vcf]
    args2=  ['-bed',outfile_bed,
             '-U', 'ALLOW_UNINDEXED_BAM']
    
    args = args1+list_bam_files+args2+call_on_this_regions
    
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command

    
def indelGenotyper_pairedSample(ref_genome, list_bam_files_normal, list_bam_files_tumor, 
                                     outfile_vcf, outfile_bed, use_mem, num_cores, path_to_gatk, 
                                     interva_list=None, call_parameters=None):
    '''
    Important Note:
    In this mode all normal/tumor bam files are joined together
    and calls is made by comparing ALL reads for normal vs tumor
    Because of that, this mode should be used with for example
    a tissue cohort, e.g Prostate, but nor with mixed cohorts,
    e.g, breast plus Prostate cohorts data. 
    '''
    # set up sample lists   
    normal_files=[]
    for bamfile in list_bam_files_normal:
        normal_files += ['-I:normal']+[bamfile]
    tumor_files=[]
    for bamfile in list_bam_files_tumor:
        tumor_files += ['-I:tumor']+[bamfile]
    
    if interva_list is not None:
        # call in specific sites
        call_on_this_regions=['-L',interva_list]
    else:
        #generate calls in all sites
        call_on_this_regions=['']

    # extra parameters
    if call_parameters is not None:
        stand_call_conf=call_parameters[0]
        stand_emit_conf=call_parameters[1]
        dcov=call_parameters[2]
    else:
        # Note 12-7-10. Check this parameters
        stand_call_conf=50.0
        stand_emit_conf=10
        dcov=50
        
    # Call
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    args1 = ['java','-Xmx'+str(use_mem)+'m', '-jar',gatk_command,'-T','IndelGenotyperV2',
            '-l', 'INFO','-R', ref_genome, '-o', outfile_vcf]
    args2=  ['-bed',outfile_bed,'-somatic',
             '-U', 'ALLOW_UNINDEXED_BAM']
    
    args = args1 + normal_files + tumor_files + args2 + call_on_this_regions
    
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command


def generate_indel_maskfile(indels_raw_bed, resources_folder, padding, indels_masked_bed):
    '''
    python python/makeIndelMask.py <raw_indels> <mask_window> <output>
    e.g.
    python python/makeIndelMask.py indels.raw.bed 10 indels.mask.bed
    '''
    mask_cmd = resources_folder+'makeIndelMask.py'
    args=['python',mask_cmd, indels_raw_bed, str(padding), indels_masked_bed]
    
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command

    
def snps_basicFiltering(ref_genome, input_raw_calls, mask_indel_file,
                       outfile_vcf, use_mem, num_cores, path_to_gatk,
                       snp_cluster_window=10):
    '''
    This function filters the snp calls according to different
    criteria:
        -quality,
        -cluster,
        -around indels
        
       java -jar /path/to/dist/GenomeAnalysisTK.jar \
      -T VariantFiltration \
      -R /seq/references/Homo_sapiens_assembly18/v0/Homo_sapiens_assembly18.fasta \
      -o /path/to/output.vcf \
      -B:variant,VCF /path/to/input.vcf \
      -B:mask,SimpleIndel /path/to/indels.calls \
      --maskName InDel \
      --clusterWindowSize 10 \
      --filterExpression "AB < 0.2 || MQ0 > 50" \
      --filterName "Nov09filters"
      -nt 6
    '''
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    args1 = ['java','-Xmx'+str(use_mem)+'m', '-jar',gatk_command,'-T','VariantFiltration',
            '-l', 'INFO','-R', ref_genome, '-o', outfile_vcf, '-B:variant,VCF',input_raw_calls]
    #basic snps filtering
    args2=  ['-B:mask,Bed',mask_indel_file,'--maskName','\"InDel\"',
             '--clusterWindowSize', str(snp_cluster_window),
             '--filterExpression', '\"MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)\"',
             '--filterName', '\"HARD_TO_VALIDATE\"'
             ]
    
    args = args1 + args2
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command


def snps_hardFiltering(ref_genome, input_raw_calls, outfile_vcf, 
                       use_mem, num_cores, path_to_gatk):
    '''
    This function filters the snp calls according to different
    criteria:
        -quality,
        -cluster,
        -around indels
        
       java -jar /path/to/dist/GenomeAnalysisTK.jar \
      -T VariantFiltration \
      -R /seq/references/Homo_sapiens_assembly18/v0/Homo_sapiens_assembly18.fasta \
      -o /path/to/output.vcf \
      -B:variant,VCF /path/to/input.vcf \
      -B:mask,SimpleIndel /path/to/indels.calls \
      --maskName InDel \
      --clusterWindowSize 10 \
      --filterExpression "AB < 0.2 || MQ0 > 50" \
      --filterName "Nov09filters"
      -nt 6
    '''
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    args1 = ['java','-Xmx'+str(use_mem)+'m', '-jar',gatk_command,'-T','VariantFiltration',
            '-l', 'INFO','-R', ref_genome, '-o', outfile_vcf, '-B:variant,VCF',input_raw_calls]

    # data specific filtering
    # Note I should allow here a better way to input non-standard filters
    # Learning from the data: covarage and other parameters
    # Consider using the AB, allelic Balance
    args3 = ['--filterExpression',
             '\"QUAL < 30.0 || DP > 40 || QD < 5.0 || HRun > 5 || SB > -0.10\"',
             '--filterName', '\"STAND_FILTER\"']
    
    
    args = args1 + args3
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command

def snps_generateVariantClusters(ref_genome,snps_calls_vcf, cluster_ouput_file,
                                 hapmap_vcf, tgk_vcf, dbsnp_vcf, path_to_gatk, use_mem):
    '''
    java -Xmx4g -jar GenomeAnalysisTK.jar \
   -R /broad/1KG/reference/human_b36_both.fasta \
   -B:input,VCF path/to/input/snpCalls.vcf \
   -B:input2,VCF path/to/another/input/snpCalls.vcf \
   -B:hapmap,VCF path/to/HapMap.vcf \
   -B:1kg,VCF path/to/1KG_project_calls.vcf 
   -B:dbsnp,VCF resources/dbsnp_132_hg18.vcf \
   -l INFO \
   -an QD -an SB -an HaplotypeScore -an HRun \
   -clusterFile path/to/output.cluster \
   -T GenerateVariantClusters
   -nt?
   '''
    # hapmap and thousand genome projects calls 
    # are taken as the true positives 
    # dbsnp is considered (by gatk developers) to have
    # lots of false positives
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    args1 = ['java','-Xmx'+str(use_mem)+'m', '-jar',gatk_command,'-T','GenerateVariantClusters',
            '-l', 'INFO','-R', ref_genome, '-clusterFile', cluster_ouput_file, '-B:input,VCF',
            snps_calls_vcf, '-B:hapmap,VCF',hapmap_vcf,'-B:1kg,VCF',tgk_vcf,
            '-B:dbsnp,VCF',dbsnp_vcf]
    
    # model covariates
    # Other covariates could be include depending on the dataset.
    args2 = ['-an', 'QD', '-an', 'SB', '-an', 'HaplotypeScore', '-an', 'HRun']

    args = args1 + args2
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command


def snps_variantRecalibration(ref_genome,snps_calls_vcf, cluster_ouput_file,
                              snps_recal_vcf, tranches_file, hapmap_vcf, 
                              tgk_vcf, dbsnp_vcf, path_to_gatk, path_to_Rscript,
                              path_to_resources,use_mem):
    '''
    java -Xmx4g -jar GenomeAnalysisTK.jar \
   -R /broad/1KG/reference/human_b36_both.fasta \
   -B:input,VCF path/to/input/snpCalls.vcf \
   -B:input2,VCF path/to/another/input/snpCalls.vcf \
   -B:hapmap,VCF path/to/HapMap.vcf \
   -B:1kg,VCF path/to/1KG_project_calls.vcf 
   -B:dbsnp,VCF resources/dbsnp_132_hg18.vcf \
   -l INFO \
   --ignore_filter HARD_TO_VALIDATE \
   --ignore_filter LowQual \
   -clusterFile path/to/output.cluster \
   -o path/to/recalibrator_output.vcf \
   -tranchesFile path/to/output.dat.tranches \
   --target_titv 2.07 \
   -resources R/ \
   -tranche 0.1 -tranche 1 -tranche 10 \
   -T VariantRecalibrator

    '''
    tmin=0.1
    tint=1
    thigh=10
    titv=3.0
    tranch_min=str(tmin)
    tranch_int=str(tint)
    tranch_high=str(thigh)
    target_titv = str(titv) # 3.0 for exomes is recommend, 2.07 for whole genomes
    
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    args1 = ['java','-Xmx'+str(use_mem)+'m', '-jar',gatk_command,'-T','VariantRecalibrator',
            '-l', 'INFO','-R', ref_genome, '-clusterFile', cluster_ouput_file, '-o',snps_recal_vcf,
            '-B:input,VCF',snps_calls_vcf, '-B:hapmap,VCF',hapmap_vcf,'-B:1kg,VCF',tgk_vcf,
            '-B:dbsnp,VCF',dbsnp_vcf]
    
    # Ignore Filters
    args2 = ['--ignore_filter', 'HARD_TO_VALIDATE', '--ignore_filter', 'LowQual']
    
    args3 =['-tranchesFile',tranches_file, '-tranche', tranch_min,'-tranche',
             tranch_int, '-tranche', tranch_high, '--target_titv', target_titv]
    args4=['-Rscript',path_to_Rscript, '-resources',path_to_resources]
    

    args = args1 + args2+ args3 + args4
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command

def snps_applyVariantCuts(ref_genome, snps_recal_vcf,snps_filtered_vcf, 
                          tranches_file, fdr_filter_level, 
                          dbsnp_vcf, path_to_gatk, use_mem):
    '''
    java -Xmx4g -jar GenomeAnalysisTK.jar \
   -R /broad/1KG/reference/human_b36_both.fasta \
   -B:input,VCF path/to/recalibrator_output.vcf \
   -B:dbsnp,VCF resources/dbsnp_132_hg18.vcf \
   -l INFO \
   --fdr_filter_level 1.0 \
   -tranchesFile path/to/output.dat.tranches \
   -o path/to/recalibrator_output.filtered.vcf \
   -T ApplyVariantCuts
    '''
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    args1 = ['java','-Xmx'+str(use_mem)+'m', '-jar',gatk_command,'-T','ApplyVariantCuts',
        '-l', 'INFO','-R', ref_genome, '-o',snps_filtered_vcf,
        '-B:input,VCF',snps_recal_vcf, '-B:dbsnp,VCF',dbsnp_vcf]
        
    args2 =['-tranchesFile',tranches_file, '--fdr_filter_level', str(fdr_filter_level)]
    

    args = args1 + args2
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command

